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ABSTRACT 


A Chebyshev-Fourier discretization with shock fitting is 
used to solve the unsteady Euler equations. The method is 
applied to shock interactions with plane waves and with a 
simple model of homogeneous isotropic turbulence. The plane 
wave solutions are compared to linear theory. 
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1. INTRODUCTION 


The effects of a shock wave upon a turbulent boundary 
layer have been studied extensively by both experimental and 
computational means. A persistent experimental result is that 
the turbulence levels downstream of the shock are higher, by 
at least a factor of 2, than the upstream levels [1]. The 
computational results, however, have consistently failed to 
exhibit this dramatic turbulence amplification [2] . This 
failure is almost surely due to inadequacies in the turbulence 
models used in the computations. Current turbulence models 
have achieved their greatest success in low-speed flows and 
attached high-speed flows vrtiich agree well with experiment 
even up to Mach 20 [3]. Evidently, some essential physics of 
separated compressible flow is missing from the models. Their 
improvement must await a better understanding of the 
interaction between shock waves and turbulence. 

Two theoretical tools are now available for the study of 
the inviscid aspects of this interaction. The basically 
linear effects may be examined analytically in a way developed 
by Ribner [4] and others, and used by Anyiwo and Bushnell [5] 
in their assessment of turbulence amplification in the shock- 
boundary layer interaction. Numerical computations based on 
the nonlinear, two-dimensional Euler equations are more 
recent. Pao and Salas [6] used this approach to study the 
shock-vortex interaction. So did Zang, Hussaini and Bushnell 
[7] in their examination of the interaction of plane waves 
with shocks. Both sets of numerical computations used second- 
order finite differences for the spatial discretization. In 
this paper we will describe a pseudospectral spatial 
discretization. 


2. NUMERICAL METHOD 

The model problem which is used to study the turbulence 
amplification and generation mechanisms is illustrated in 
Figure 1. At time t = 0 an infinite, normal shock at x = 0 
separates a rapidly moving, uniform fluid on the left from the 
fluid on the right vrtiich is in a quiescent state except for 
some specified fluctuation. The initial conditions are chosen 
so that in the absence of any fluctuation the shock moves 
uniformly in the positive x-direction with a Mach number 
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(relative to the fluid on the right) denoted by M . i n the 
presence of fluctuations the shock front will develop 
ripples. The structure of the shock is described by the 
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Figure 1: Model problem in the physical domain 
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function x g (y,t). The numerical calculations are used to 
determine the state of the fluid in the region between the 
shock front and some suitable left boundary, x^(t), and also 
to determine the motion and shape of the shock front itself . 

In the applications given below the physical problem is 
periodic in y with period y L . The change of variables 


x - x T (t) 
X o(y»t) - Xr (t) 


Y = y/y L 

T = t 

produces the computational domain 

0 < X < 1 
0 < Y < 1 
T > 0. 


( 1 ) 


( 2 ) 


In terms of the computational coordinates the two-dimensional 
Euler equations are 


Q t + BQ X + CQ y = 0, 


where Q = (P,u,v,S), 
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(5) 


The contravariant velocity components are given by 

U = X t + uX x + vXy 


and 


V = Y t + uY x + vY y . 


( 6 ) 


A subscript denotes partial differentiation with respect to 
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the indicated variable. P, c and S are the natural 
logarithm of pressure, the sound speed and the entropy 
(divided by the specific heat at constant volume), 
respectively, all normalized by reference conditions at 
downstream infinity; u and v are velocity components in 

the x and y directions, both scaled by the characteristic 
velocity defined by the square root of the pressure-density 
ratio at downstream infinity. The ratio of specific heats is 
denoted by Y; a value Y = 1.4 has been used for all the 
calculations in this paper. 

Let n denote the time level and At the time 

increment. The time discretization of Eq. (3) is 

Q = [l-AtL n ]Q n (7) 

Q n+1 = \ [Q n +(l-AtL)Q], (8) 

where L denotes the spatial discretization of B3 X + CSy* 
The solution Q has the Chebyshev - Fourier series expansion 

M N/2-1 

Q(X,Y,T) = l l Q <T)t (£)e 2iriqY , (9) 

p=0 q=-N / 2 pq p 


where 

degree 


5 = 2X-1 and T m is the Chebyshev polynomial of 
m. The derivatives Q x and Qy are approximated by 


M 



T‘ 1 q (1 - 0) (t)t 

q=-N/2 Pq P 


( 10 ) 


where 


with 


and 


M 

Qy - 2TT l 

p=0 


q— N/2 Pq P 


( 11 ) 
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pq 
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( 12 ) 


(13) 


Q 


( 0 , 1 ) 

pq 


iq w 


(14) 


The most critical part of the calculation is the 
treatment of the shock front. The shock-fitting approach used 
here is desirable because it avoids the severe post-shock 
oscillations that plague shock-capturing methods. The time 
derivative of the Rankine-Hugoniot relations provides an 
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equation for the shock acceleration. This equation is 
integrated to update the shock position. The jump conditions 
themselves are then used at x g (y,t) to determine the 
variables on the left (high perssure) side of the shock from 
the specified downstream flow (see [8] for details). For 
M g > 2.08, the left boundary at X L (t) is a supersonic inflow 
boundary. Lower Mach numbers require more sophisticated 
boundary conditions as discussed in [8]. Spectral methods for 
compressible flows typically require some sort of filtering. 


3. INTERACTION OF PLANE WAVES WITH A SHOCK 

The non-linear interaction of plane waves with shocks was 
examined at length in [7] . The numerical method used there 
was similar to the one described above but employed second- 
order finite differences in place of the present 
pseudospectral discretization. Detailed comparisons were made 
in [7] with the predictions of linear theory [4] . The linear 
results turned out to be surprisingly robust, remaining valid 
at very low (but still supersonic) Mach numbers and at very 
high incident wave amplitudes. The only substantial 
disagreement occurred for incident waves whose wave fronts 
were nearly perpendicular to the shock front. This type of 
shock-turbulence interaction is a useful application of the 
pseudospectral technique. The method can be calibrated in the 
regions for which linear theory has been shown to be valid. 
The questionable regions can then be gainfully investigated by 
means of this alternative method. 

We concentrate on Mach 3 normal shocks interacting with 
incident vorticity waves. The perturbation velocities ahead 
of the shock are taken to be 

(u',vO = A"(-k 1>y ,k 1>x )(c 1 /p 1 k 1 )cos(k 1 »j t ), 

where the wavevector «* (k^ x ,kj ) and k' is the 
amplitude. In terms of the incidence angle 6j, = 

(kjcos 0j,k^sin 9jJ. The transmitted vorticity behind the 
shock can be expressed in the same manner with all 
subscripts 1 changed to 2. The vorticity transmission 
coefficient is defined to be the ratio A^/A^. At Mach 3 the 
ratio of the rotational part of the perturbation velocities 
is 0.158 times the transmission coefficient; if these 
perturbation velocities are scaled by the respective mean 
velocities (in the frame in which the shock is at rest) , then 
this ratio is 0.611 times the transmission coefficient. 

A typical numerical simulation is shown in Figure 2. At 
time t = 0 all the fluid in the computational domain, which 
then extended from x^ = -0.5 to x g = 0.0, was in uniform 




Figure 2: Responses at time t = 0.80 to a 0.1% 
vorticity wave striking a Mach 3 shock at a 
30° angle of incidence. The shock front is 

denoted by the solid line connecting the top 
and bottom boundaries. Dashed lines indicate 
negative contour levels. 



7 


flow. The flow ahead of the shock was given by Eq. (15) with 
an amplitude k' = 0.001. In order to reduce the transients 
that would result If the fluid were to encounter this wave in 
a sudden fashion, the amplitude was turned on smoothly from 
k' = 0.000 to A" = 0.001 during the passage of the shock 
over the first half -wavelength of the incident wave. The 
interaction of the shock with the incident vorticity wave has 
produced a transmitted vorticity wave as well as generated 
acoustic and entropy waves. The acoustic wave clearly has a 
larger Inclination angle and a faster propagation speed than 
the other two waves. 



X 


Figure 3: Post-shock dependence at time t « 0.80 of the 

vorticity response to a vorticity wave incident 
at 30° to a Mach 3 shock. The solid line is 
the linear theory prediction. Computed 

responses are given for 0.1% (circles) and 
50% (diamonds) waves. 

Figure 3 indicates how we measure the transmitted 
vorticity. The computed velocities are differentiated 
pseudospectrally and combined to produce the response 
vorticity at each point on the computational grid. At each 
fixed value of X we perform a Fourier analysis in Y of the 
vorticity. . The wavenumber of interest is q = 1 (see Eq . 
(9)). This amplitude is then scaled as indicated above to 

yield A£ and thence the ratio A^/A^ which is plotted in 
Figure 3 as a function of a typical value of the physical 

coordinate x corresponding to X. In the case shown in 

Figure 3 the incident wave does not reach full strength 
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until t = 0.56. A single value representing the transmission 
coefficient is obtained by averaging the x-dependent 
responses. The standard deviation of these values serves as 
an error estimate, indicated by the error bars in the 
following two figures. 

The angular dependence of the vorticity transmission 
coefficient for Mach 3 is displayed in Figure 4. Correspond- 
ing results from the finite difference discretization used in 
[7] are shown for comparison. The amplitude dependence for 
30° incident waves is shown in Figure 5. 



Figure 4: Incident angle dependence of the vorticity 
response to 0.1% vorticity waves incident 
upon a Mach 3 shock. Pseudospectral (circles) 
and finite difference (squares) results are 
given. The solid line is the linear theory 
prediction for sub-critical angles. 

Linear theory predicts that the vorticity responses are 
peaked in the vicinity of 0^ = 62°, the so-called critical 
angle. At larger angles the acoustic responses are 
evanescent. These predictions do not appear in the non-linear 
calculations. Tentative explanations for this failure of 
linear theory are provided in [7]. 

Below the critical angle, however, the linear theory is 
extremely robust. Figure 5 indicates that it remains valid 
for incident relative velocity fluctuations aproaching 50%. 
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(In the frame in which the shock is at rest, the ratio of the 
fluctuating velocity to the mean velocity is one-third of the 
amplitude • ) 



Figure 5 Amplitude dependence of the vorticity response 
to vorticity waves incident at 30^ to a Mach 
3 shock. The solid line is the linear theory 
prediction. 


4. INTERACTION OF TURBULENCE WITH A SHOCK 

Low intensity turbulence can be viewed as the 
superposition of plane shear waves just discussed. Figure 6 
shows velocity vector plots of the interaction of a shock with 
a simple homogenous isotropic turbulence model in which the 
velocities are defined as sums of sinusoidal waves with random 
phases. The frame at t = 0.38 shows the undisturbed 
velocity field downstream of the shock and t = 1.6 shows the 
result after the shock has traveled approximately one 
wavelength of the longest wave component of the undisturbed 
turbulence. The development of the post-shock vortices seen 
in the figure coincides with the generation of rows of 
pressure and entropy spots behind the shock. 


5. CONCLUSION 

We have solved the two-dimensional unsteady non-linear 
Euler equations with a Chebyshev-Fourier pseudospectral 
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method. The linear theory of the Interaction of a shock and 
plane shear waves has provided a test of the method. Details 
of the smoothing methods and questions of accuracy will appear 
in a later paper. But accuracy comparable to the finite 
difference calculations can be obtained with about half the 
number of Chebyshev modes in X and with only four modes in 
Y. The results indicate that it should be possible to compute 
more complicated shock-turbulence interactions. 



Figure 6: Velocity vector plots at t = 0.38 (top) and 

t = 1.6 (bottom) of the interaction of a Mach 
3 shock with low intensity homogeneous 
isotropic turbulence. 
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